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Abstract. The efficient numerical soiution of Non-LTE muitiievei transfer problems 
requires the combination of highly convergent iterative schemes with fast and accurate 
formal solution methods of the radiative transfer (RT) equation. This contribution^ be- 
gins presenting a method for the formal solution of the RT equation in three-dimensional 
(3D) media with horizontal periodic boundary conditions. This formal solver is suitable 
for both, unpolarized and polarized 3D radiative transfer and it can be easily combined 
with the iterative schemes for solving non-LTE multilevel transfer problems that we have 
developed over the last few years. We demonstrate this by showing some schematic 3D 
multilevel calculations that illustrate the physical effects of horizontal radiative trans- 
fer. These Non-LTE calculations have been carried out with our code MUGA 3D, a 3D 
multilevel Non-LTE code based on the Gauss-Seidel iterative scheme that Trujillo Bueno 
and Fabiani Bendicho (1995) developed for RT applications. 
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1. Introduction 

To what extent can we trust diagnostic results obtained with the assump- 
tion that the solar atmospheric plasma is composed of homogeneous plane- 
parallel layers or via approximations that neglect horizontal radiative trans- 
fer (RT) effects? How important are the errors in the magnetic fields, tem- 
peratures and velocities inferred by confronting spectro-polarimetric ob- 
servations with Non-LTE ID RT model calculations? Clearly, to provide 
proper answers to questions like these requires to develop first efficient 3D 
RT methods that allow Non-LTE effects in complex atomic models with 
many levels and transitions to be rigorously investigated. 

^Published in 1999 in the book Solar Polarization, edited by K.N. Nagendra & J.O. 
Stenflo. Kluwer Academic Publishers, 1999. (Astrophysics and Space Science Library ; 
Vol. 243), p. 219-230 
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There is a second reason which makes the development of fast iterative 
methods for 3D Non-LTE RT so relevant. This is because processes of 
energy exchange by radiation play an important role in the structure and 
dynamical behaviour of the stellar magnetized plasma. Thus, for instance, if 
one wishes to perform time-dependent radiation hydrodynamics simulations 
similar to those carried out by Carlsson and Stein (1997), but in 3D instead 
of ID, it turns out to be imperative to have first access to numerical methods 
capable of accurately yielding the self-consistent atomic level populations 
at the cost of only very few formal solution times. 

The efficient solution of multilevel transfer problems requires the com- 
bination of a highly convergent iterative scheme with a fast formal solver of 
the RT equation. In Section 2 we briefly comment on a hierarchy of iterative 
schemes that can be applied for solving multilevel Non-LTE problems with 
increasing improvements in the convergence rate and total computational 
work. The 3D multilevel transfer calculations that we present in this con- 
tribution have been obtained by combining a highly convergent iterative 
scheme based on Gauss-Seidel iteration (Trujillo Bueno and Fabiani Ben- 
dicho, 1995) with a fast 3D formal solver that has parabolic accuracy (see 
Section 3). As was the case with our 2D formal solver, our generalization 
to 3D is based on the "short-characteristics" method of Kunasz and Auer 
(1988). Our 3D multilevel code is called MUGA 3D ("Multi-level Gauss- 
Seidel Method") and it is substantially faster than our code MALI 3D, 
which is based on Jacobi iteration (see Rybicki and Hummer, 1991; Auer, 
Fabiani Bendicho and Trujillo Bueno, 1994). 

Section 3 briefly describes our 3D formal solver as applied to the scalar 
transfer equation for the specific intensity (I). In order to be able to consider 
3D atmospheric models where solar plasma structures repeat themselves 
along the horizontal directions we choose horizontal periodic boundary con- 
ditions along the Cartesian coordinates X and Y. Although we do not give 
any details here, we have also generalized to 3D the Stokes-vector ID for- 
mal solver method developed by Trujillo Bueno (1998), which is based on 
the matrix exponential approximation to the evolution operator. 

In Section 4 we show some illustrative 3D multilevel transfer calculations 
for a 5-level Ca II model atom where the H, K and infrared triplet lines are 
treated simultaneously, taking fully into account the interlockings by which 
photons are converted back and forth between the different line transitions 
in the assumed 3D medium. Here we consider schematic 3D solar models 
characterized by horizontal sinusoidal temperature inhomogeneities. With 
the help of these 3D multilevel calculations we are able to illustrate some 
subtle effects of horizontal radiative transfer that are important for the cor- 
rect interpretation of high spatial resolution observations. Finally, Section 
5 gives our conclusions. 



THREE-DIMENSIONAL RADIATIVE TRANSFER 



3 



2. Iterative Methods for Multilevel Transfer 

The simplest procedure one might think of to solve self-consistently the 
kinetic and RT equations is A— iteration: using the current estimate of the 
atomic level populations at each spatial grid-point (or, more in general, of 
the irreducible tensor components of the atomic density matrix; see Trujillo 
Bueno 1999) evaluate the absorption and emission coefficients. Next, solve 
the radiative transfer equation and compute the radiation field intensity 
in all transitions. Finally, with the radiative rates obtained, solve the ki- 
netic equations at each spatial grid-point independently and obtain a new 
estimate of the atomic level populations. However, as is well known, under 
typical NLTE conditions in optically thick media this A— iteration method 
converges extremely slowly. This is certainly regrettable because with this 
method there is no need of inverting large matrices and the computing time 
per iteration is minimal. In any case, as demonstrated by Trujillo Bueno and 
Manso Sainz (1999), in solar-like atmospheres the A— iteration method can 
be used to find the self-consistent solution of Non-LTE polarization transfer 
problems if one initializes using the "exact" solution corresponding to the 
unpolarized transfer case. 

The dream of numerical RT is to develop iterative methods where ev- 
erything goes as simply as with the A— iteration scheme, but for which the 
convergence rate is extremely high. The iterative methods for RT appli- 
cations based on Gauss-Seidel iteration that we have developed over the 
last few years have been worked out with this motivation in mind. Their 
convergence rate is extremely high, there is no need of constructing and 
inverting any large matrix, and the computing time per iteration is similar 
to that of the A— iteration method. A full account of these developments 
can be found in the following publications: 

1) Auer, Fabiani Bendicho and Trujillo Bueno (1994) present the gener- 
alization of the Jacobi-based MALI method of Rybicki &: Hummer (1991) 
to multilevel RT in 2D. They also developed a short-characteristics strategy 
to do the formal solution in 2D Cartesian coordinates with horizontal pe- 
riodic boundary conditions. Of particular interest is a simple grid- doubling 
technique which both rapidly finds the converged solution in fine meshes 
and automatically estimates its corresponding true error. The total compu- 
tational work scales as NP^ , with NP the total number of spatial grid-points 
in a fixed computational domain. 

2) Trujillo Bueno and Fabiani Bendicho (1995) developed a novel iter- 
ative scheme based on Gauss-Seidel (GS) iteration (MUG A). This is the 
paper on which our present work is based on. The total computational work 
scales as NP^/4 for the pure GS method (implemented as suggested in the 
conclusions of their paper), and as NPy^NP for the successive overrelaxation 
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(SOR) method. This paper was fundamental for a successful development 
of our non-linear multigrid method for RT applications. 

3) Fabiani Bendicho, Trujillo Bueno and Auer (1997) consider the ap- 
plication of the non-linear multigrid method (see Hackbush, 1985) to mul- 
tilevel RT. Here the iterative scheme is composed of two parts: a smoothing 
one where a small number of MUGA iterations on the desired finest grid 
are used to get rid of the high-frequency spatial components of the error 
in the current estimate, and a correction obtained from the solution of an 
error equation in a coarser grid. With this method the total computational 
work scales simply as NP, although it must be said that the computing time 
per iteration is about 4 times larger than that required by the A— iteration 
method. 
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Figure 1. Variation with the grid-spacing of the maximum eigenvalue of the iteration 
operator corresponding to several multilevel iterative schemes. The MG symbol refers to 
our non-linear multigrid code. 

In order to compare the convergence rates of these three iterative meth- 
ods, we present in Fig. 1 an estimate of the maximum eigenvalues (p) of 
the corresponding iteration operator, which controls the convergence prop- 
erties of such iterative schemes (see Trujillo Bueno and Fabiani Bendicho, 
1995). The knowledge of this maximum eigenvalue (p) is useful because 
errors decrease as p***", where "itr" is the iterative step. We obtain this 
information from multilevel Ca II calculations in several ID grids with de- 
creasing grid-size (Az) by calculating Rc{itr -\- 1)/Rc{itr) for itr » 1, 
where Rc{itr) is the maximum relative change in the level populations. As 
it can be noted in Fig. 1 the convergence rate of both, the MALI and MUGA 
schemes decreases when the spatial resolution of the grid is improved, while 
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Figure 2. 3D Cartesian spatial grid surrounding the grid-point of interest (O) where the 
specific intensity lo is to be calculated. This is done by solving analytically the integral 
of Eq. (2) along the short- characteristics MO corresponding to the ray of direction f2. 



the maximum eigenvalue of our non-linear multigrid method is always very 
small (p ~ 0.1) and insensiUve to the grid-size. A maximum eigenvalue 
p = 0.1 means that the error decreases by one order of magnitude each 
time we perform an iteration. This explains that, typically, two multigrid 
iterations are sufficient to reach the self-consistent solution for the atomic 
level populations. 

The 3D multilevel calculations shown in this contribution were obtained 
with our code MUGA 3D, i.e. with a multilevel GS scheme based on the 
paper by Trujillo Bueno and Fabiani Bendicho (1995) combined with the 
following 3D formal solver. In practical applications we always use MUGA 
3D with Ng (1974) acceleration. 



3. The 3D formal solver. 

The scalar RT equation for the specific intensity is 

^=X.(S. -I.), (1) 

where s is the geometric distance along the ray propagating in a certain 
direction in a 3D medium, Xv is the total opacity and Sj^ the source function. 

We now consider a 3D Cartesian spatial grid (see Fig. 2). Point O is 
the grid-point of interest at which one wishes to calculate the specific in- 
tensity To, for a given frequency {v) and a direction (p). Point M is the the 
intersection point with the grid-plane that one finds when moving along 
— O. At this upwind point we assume that the specific intensity Im for the 
same frequency and angle is known from previous steps. In a similar way, 
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point P is the intersection point with the grid-plane that one encounters 
when moving along Q. We also introduce the optical depths along the ray 
between points M and O (Atm) and between points O and P (Arp). Prom 
the formal solution of the previous transfer equation one finds that 

Io = lMe-^^+/ S(t)e-(^^-*)dt, (2) 
Jo 

where the optical depth variable t is measured from M to O. 

The integral of Eq. (2) can be solved analytically by integrating along 
the short- characteristics MO assuming some prescribed variations for the 
source function (e.g. linear variation along M and O, or parabolic along 
M,0 and P, or cubic-ccntcrcd around point O, etc.). Our 3D formal solution 
method assumes that the source function S{t) varies parabolically along M,0 
and P. The result reads: 

lo = iMe-^^ + *mSm + *oSo + ^'pSp, (3) 

where ^'x (with X either M, O or P) are given in terms of the quantities 
Atm and Arp that we evaluate numerically by assuming that ln(x) varies 
linearly with the geometrical depth, x being the opacity. 

It is very important to point out that one should avoid the use of a 
formal solution method based on a linear interpolation formula, i.e. one 
should avoid assuming that, for each grid-point O of interest, the source 
function varies linearly along points M and O. Otherwise, the accuracy of 
the self-consistent solution will never be better than about 10%, even by 
choosing a very large number of grid-points per opacity scale height (see 
Trujillo Bueno, 1998). The reason for this is that the use of linear inter- 
polation for S{t) leads to formal solution methods that are unable to yield 
the correct asymptotic behaviour for the intensity when having nonlinear 
source functions in optically thick atmospheres. We emphasize that our 3D 
formal solution method is based on the above-mentioned parabolic approx- 
imation and it only uses the linear approximation formula for calculating 
the radiation field at the upper and lower boundaries for rays going out of 
such boundaries. However, as we illustrate below, the use of the parabolic 
approximation for investigating problems where we have sudden variations 
in the physical quantities requires to implement it using an improved ver- 
sion of the monotonic upwind interpolation technique applied by Auer and 
Paletou (1994). 

The application of this formal solution method in ID is straightforward. 
A detailed description of how to implement it in 2D slabs with prescribed 
irradiation on the lateral boundaries can be found in Kunasz and Auer 
(1988) and Auer and Paletou (1994). In 2D with horizontal periodic bound- 
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ary conditions is slightly more complicated and a suitable strategy has been 
described by Auer, Fabiani Bendicho and Trujillo Bueno (1994). 

The main changes when going to 3D imposing horizontal periodic bound- 
ary conditions lie in the interpolation. We have assumed that Im is known 
(see Fig. 2) but, in most cases, the M-point (hke the point P) will not be a 
grid-point of the chosen 3D spatial grid. The intensity at this M-point has 
to be calculated by interpolating from the available information at the nine 
surrounding grid-points, as we must also do for obtaining the opacities and 
source functions at M and P. 

Parabolic interpolation can however generate spurious negative inten- 
sities if the spatial variation of the physical quantities is not well resolved 
by the spatial grid. This happens, for instance, if one tries to simulate 
the propagation of a beam in vacuum using a three dimensional grid. The 
analytical solution in this case is simply Iq = Im) and any numerical er- 
ror is due to the effect of the conventional way of applying the parabolic 
interpolation. To avoid these problems we have improved the ID mono- 
tonic interpolation strategy of Auer and Paletou (1994), and generalized it 
to the two-dimensional parabolic interpolation that is required for 3D RT 
calculations (Fabiani Bendicho and Trujillo Bueno, in preparation). 




Figure 3. Ray propagating in vacuum. Emergent intensity with parabolic interpolation 
implemented without (a) and with (b) our improved version of the monotonic upwind 
interpolation strategy. 



In Fig. 3 we show the intensity that emerges from a 3D computational 
box after having illuminated its lower boundary with a beam of given in- 
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tensity. Here the aim is to simulate a beam propagating in vacuum. Using 
standard parabolic interpolation (Fig. 3a) leads to unphysical negative in- 
tensities. However, with our improved version of the upwind monotonic 
interpolation strategy we guarantee that the parabolic interpolation is not 
introducing spurious sources and sinks and its second order accuracy is 
maintained. 

For solving Non-LTE radiative transfer problems as shown above the 
chosen formal solution method has to be implemented in a way such that 
it rapidly computes, at each spatial grid-point, the radiation field for each 
frequency and direction of the chosen numerical quadratures together with 
the diagonal element of the A— operator. To this end, one may also try 
with higher-order methods if desired. We have found that our parabolic 
formal solvers are fast and accurate enough, stable and easy to work with 
independently of the geometry and of the iterative method used. 

4. Illustrative 3D multilevel results 

Our aim here is to study a few illustrative results for Ca H in a schematic 3D 
atmospheric model characterized by the following temperature structure: 

T(X,Y,Z) = 5000 + 500sin(aX) sin(bY), (4) 

where a = 27r/Px and b = 27r/Py, with Pj; and Py the horizontal periods 
along the X and Y Cartesian coordinates, respectively. For this example 
we choose Px = Py = 1000 km. We also assume that the total hydrogen 
number density is exponentially stratified along the vertical direction (Z) 
with a scale height 7^=100 km. We adopt a standard 5-level Ca H atomic 
model and a constant electron density rig = 10^^ cm~^. We compare our 
full 3D multilevel results with the corresponding 2D multilevel calculation 
and with results obtained using a well-known approximation that neglects 
the effect of horizontal RT on the level populations, i.e. with the so-called 
1.5D approximation (Mihalas, Auer and Mihalas, 1978; Solanki, Steiner 
and Uitenbroek, 1991). 

We point out that our calculations do take into account the change in 
the line opacities caused by the assumed 500 K horizontal temperature in- 
homogcneities. Thus, the 3D multilevel transfer effects are the result of the 
combined action of the sm,oothing effect of the source-function fluctuations 
(Spiegel, 1957; Kneer and Trujillo Bueno, 1987) and the channelling effect 
due to the opacity fluctuations (Cannon, 1970; Trujillo Bueno and Kneer, 
1990). Since we use a 5-level atomic model, we can examine the 3D effects 
for the H, K and the infrared triplet lines. Figures 4, 5 and 6 give the vari- 
ation with X and Y of the vertically emergent intensity at the continuum 
frequency (upper part) and at the line core (lower part). We only show this 
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for the H and 8662 A lines since for the K hne and the other two infrared 
lines very similar results were found. 




Figure 4- The vertically emergent intensity at the line core (lower surface plot) and at 
the continuum frequency (higher surface plot) corresponding to the 1.5D approximation 
applied to a model atmosphere with = Py = 1000 km. These intensities are shown 
for the H and 8662 A Ca II lines in the left and right panels of the figure, respectively. 
Geometrical distances are measured in units of the opacity scale height Ti. 

Fig. 4 corresponds to calculations carried out with the 1.5D approxi- 
mation. As pointed out above this approximation completely neglects the 
effects of horizontal RT on the level populations. Consequently, the result 
is that, not only at the continuum frequency, but also at the hne cen- 
ter the emergent intensity shows large horizontal fluctuations. The 1.5D 
atomic level populations have large horizontal fluctuations at all atmo- 
spheric depths, simply because one is here neglecting the smoothing and 
channelling horizontal 3D RT effects when solving the kinetic and RT equa- 
tions. 

Fig. 5 shows the 2D multilevel transfer results (see also Auer, Fabi- 
ani Bendicho and Trujillo Bueno, 1994). Here it was assumed that the 
model's temperature was only fluctuating along the horizontal X-direction 
with Px = 1000 km; i.e. horizontal transfer effects along the Y-direction 
are completely neglected. The explanation of the results is the following: 
above the height of thermalization, 2D horizontal transfer effects are effi- 
ciently damping the horizontal fluctuations of the populations of the two 
uppermost levels of the chosen 5-level Ca II model atom. This smoothing 
is translated to the line source functions, which are proportional to n^A^^i 
for the H line and to 774^4^2 for the 8662 A line, with A the Einstein coef- 
ficient for spontaneous emission. Thus, since the continuum photons come 
from the deepest layers, one sees strong horizontal intensity fluctuations 
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Figure 5. 2D vertically emergent intensities. For more information see Fig. 4. 

at continuum frequencies, while at the line core one finds horizontal inten- 
sity fluctuations of very small amplitude because these photons come from 
the higher atmospheric layers. Note that, both for the H and the 8662 A 
infrared lines, the vertically emergent intensities are fluctuating in phase 
with the assumed thermal inhomogeneities. 




Figure 6. 3D vertically emergent intensities. For more information see Fig. 4. 

Fig. 6 depicts the results of our full 3D multilevel calculation with Px = 
Py = 1000 km. For the H line the result is qualitatively similar to the 2D 
results, although the amplitude of the horizontal intensity fluctuations at 
line-center is significantly smaller, as expected from the fact that in 3D the 
horizontal transfer effects are enhanced with respect to those found for the 
2D case. However, for the 8662 A line we find a result that, at first sight. 



THREE-DIMENSIONAL RADIATIVE TRANSFER 



11 



might appear surprising: at line center the resulting horizontal fluctuation 
of the vertically emergent intensities is anticorrelated with respect to the 
fluctuation of the assumed thermal inhomogeneities. In other words, at 
the line core of this infrared line, one gets more intensity from the coolest 
atmospheric points than from the hottest ones. The explanation is that the 
IR line source function shows a change, by vr, in phase above the height (he) 
in the atmosphere where the amplitude of the horizontal fluctuation in 714 is 
greatly diminished. And this change of phase is, in turn, due to the fact that 
the line source function is the result of dividing an almost non-fluctuating 
emissivity by a line opacity (essentially set by 71,2) that fluctuates in phase 
with the assumed temperature inhomogeneities. This occurs both in 2D 
and in 3D, but with the difference that h'^^ > h^^ for given values of Px 
and Py. In fact, one finds indeed a similar result in 2D, but for Px values 
significantly smaller than 1000 km. Since the horizontal transfer effects are 
more important in 3D than in 2D, it turns out that such an effect already 
sets in at Px = Py = 1000 km. 

We thus see that our 3D multilevel calculations demonstrate the effi- 
ciency of horizontal RT effects and confirm (sec Kneer, 1981; Trujillo Bueno, 
1990) that the interpretation of high-spatial resolution observations ignor- 
ing the existence of horizontal transfer effects (as is done when using the 
1.5D approximation) may substantially underestimate the actual inhomo- 
heneities present in the solar atmospheric plasma. 

5. Concluding Remarks 

We have developed a 3D multilevel code (MUGA-3D) that combines the 
Gauss-Seidel iterative scheme of Trujillo Bueno and Fabiani Bendicho (1995) 
with a 3D formal solver that uses horizontal periodic boundary conditions. 
With this new code we have performed some 3D multilevel simulations that 
highlight the importance of carefully investigating the effects of horizontal 
radiative transfer using realistic atmospheric and atomic models. 

We point out that our 3D formal solver can be used not only for solving 
"unpolarized" multilevel transfer problems, but also resonance line polar- 
ization and Hanle effect problems, like those considered in these Proceed- 
ings by Manso Sainz and Trujillo Bueno (1999), Paletou et. al. (1999) or 
Dittmann (1999). This is because, for these polarization transfer cases, the 
absorption matrix is diagonal. As a result, we have similar equations for 
the Stokes I, Q and U parameters. 

However, for the solution of more general polarization transfer problems, 
like the Non-LTE Zeeman line transfer case considered by Trujillo Bueno 
and Landi Degl'Innocenti (1996), but in 3D instead of ID, one needs a 
3D formal solution method of the Stokes-vector transfer equation. This is 
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because here the absorption matrix turns out to be a full 4x4 matrix and all 
the Stokes parameters are coupled together. To this end we have generalized 
to 3D the Stokcs-vcctor formal solver developed by Trujillo Bueno (1998), 
which can be considered as a generalization to polarization transfer of the 
short-characteristics method. 

We would like to end this paper by saying that over the last 10 years we 
have witnessed impressive progress concerning the development of highly 
convergent iterative schemes and accurate formal solvers for RT applica- 
tions. Now it is time to apply them with physical intuition in order to 
improve our knowledge of the Sun, its magnetic field and its polarized 
spectrum. 
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